Beta diversity

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("openxlsx", "kableExtra", "dynamicTreeCut", "cowplot", "ggpubr")
invisible(lapply(packages, require, character.only = TRUE))

HCA on raw data

Preparation

# load rdata file with phyloseq object
setwd(path_rdata)
ps <- readRDS("MED_phyloseq.rds")

# transform count table into percent table
ps.percent <- transform_sample_counts(ps, (function(x) x / sum(x)))

# extract otu and tax tables
otu <- data.frame(ps.percent@otu_table)
tax <- data.frame(ps.percent@tax_table)
sam <- data.frame(ps.percent@sam_data)

# add Genus column in otu table
otu$Genus <- tax$Genus
otu <- otu %>% select(c(Genus, everything()))

# select Wolbachia
otu_wolbachia <- otu[otu$Genus=="Wolbachia" & !is.na(otu$Genus),] %>% t() 
otu_wolbachia <- otu_wolbachia[-1,]

# add sample names as rownames
samples_names <- rownames(otu_wolbachia)

# convert column into numeric type
otu_wolbachia <- data.frame(apply(otu_wolbachia, 2, function(x) as.numeric(as.character(x))))

# add column Wolbachia with sum of Wolbachia percentage
otu_wolbachia$Wolbachia <- rowSums(otu_wolbachia)

# set rownames 
rownames(otu_wolbachia) <- samples_names

# add column Sample with sample id
otu_wolbachia$Sample <- rownames(otu_wolbachia)
otu_wolbachia <- otu_wolbachia %>% select(c(Sample, Wolbachia))

# print
otu_wolbachia %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Wolbachia
CTC1 CTC1 0.0000000
CTC10 CTC10 0.0352626
CTC11 CTC11 0.0077123
CTC12 CTC12 0.0191972
CTC13 CTC13 0.0026615
CTC14 CTC14 0.1074139
CTC15 CTC15 0.0094047
CTC2 CTC2 0.0087316
CTC3 CTC3 0.0007014
CTC4 CTC4 0.0088827
CTC5 CTC5 0.0017734
CTC6 CTC6 0.0020249
CTC7 CTC7 0.0050378
CTC9 CTC9 0.0058233
NP14 NP14 0.9888611
NP2 NP2 0.9993923
NP20 NP20 0.4705882
NP27 NP27 0.9789303
NP29 NP29 0.0492611
NP30 NP30 0.2894737
NP34 NP34 0.3157895
NP35 NP35 0.0032998
NP36 NP36 0.0000000
NP37 NP37 0.0002480
NP38 NP38 0.0218636
NP39 NP39 0.0023176
NP41 NP41 0.0000162
NP42 NP42 0.0004591
NP43 NP43 0.0070779
NP44 NP44 0.0000127
NP5 NP5 0.9999973
NP8 NP8 0.9999492
S100 S100 1.0000000
S101 S101 0.8433759
S102 S102 1.0000000
S104 S104 0.9999237
S105 S105 0.9999100
S106 S106 1.0000000
S107 S107 0.9999808
S108 S108 0.9999641
S109 S109 0.9999492
S110 S110 1.0000000
S113 S113 0.9941796
S121 S121 1.0000000
S122 S122 1.0000000
S123 S123 1.0000000
S124 S124 0.9999449
S126 S126 0.9997374
S127 S127 1.0000000
S128 S128 1.0000000
S146 S146 0.9993204
S147 S147 0.9999205
S148 S148 0.9993646
S150 S150 1.0000000
S151 S151 1.0000000
S152 S152 0.9998674
S153 S153 0.9998826
S154 S154 0.9993767
S160 S160 0.9999035
S162 S162 0.9999603
S163 S163 0.9999188
S164 S164 0.9993741
S165 S165 0.9995488
S166 S166 0.9992847
S167 S167 0.9997153
S169 S169 1.0000000
S170 S170 1.0000000
S18 S18 0.9909091
S19 S19 0.6928881
S20 S20 0.9376633
S21 S21 0.6642608
S22 S22 0.5062214
S23 S23 0.9421566
S24 S24 0.9054080
S25 S25 0.9970647
S26 S26 0.9733037
S27 S27 0.2155884
S28 S28 0.9646902
S29 S29 0.0098804
S30 S30 0.1662131
S31 S31 0.1847797
S32 S32 0.2115416
S33 S33 0.3796025
S34 S34 0.5882353
S35 S35 0.4507416
S36 S36 0.1426408
S37 S37 0.0420938
S38 S38 0.8371053
S39 S39 0.3274431
S40 S40 0.4336303
S41 S41 0.3043478
S42 S42 0.9997565
S43 S43 0.9369888
S44 S44 0.9795664
S45 S45 0.9815447
S46 S46 0.4921260
S47 S47 0.4081812
S48 S48 0.9938841
S49 S49 0.9925667
S50 S50 0.9898880
S51 S51 0.9936881
S52 S52 0.9833898
S53 S53 1.0000000
S55 S55 0.9999604
S56 S56 1.0000000
S57 S57 0.9998373
S58 S58 0.9977978
S59 S59 1.0000000
S60 S60 1.0000000
S61 S61 1.0000000
S63 S63 0.9955108
S64 S64 1.0000000
S65 S65 0.3235294
S77 S77 0.0000000
S79 S79 1.0000000
S80 S80 1.0000000
S83 S83 1.0000000
S84 S84 1.0000000
S85 S85 0.9994486
S86 S86 1.0000000
S87 S87 0.9998787
S88 S88 0.9977787
S89 S89 0.0032180

HCA

Run HCA

# Compute distance on samples with Wolbachia
d.wolbachia <- dist(otu_wolbachia[,-1, drop=FALSE], method="euclidean")

# HCA
hca.ward <- hclust(d.wolbachia,method="ward.D2")

Dendogram

# 2 groups
plot(hca.ward)
rect.hclust(hca.ward,k=2)

# 3 groups
plot(hca.ward, ask=FALSE)
rect.hclust(hca.ward, k=3)

# set displaying
rh <- rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))

beg_clus <- head(cumsum(c(1, lengths(rh))), -1)
beg_clus <- beg_clus + 5
y_clus <- weighted.mean(rev(hca.ward$height)[2:3], c(4, 1))

# plot
plot(hca.ward, cex=0.7) 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)

# save plot
setwd(path_plot)
tiff("1Ec_dendogram.tiff", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="") 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()
## quartz_off_screen 
##                 2
png("1Ec_dendogram.png", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="") 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()
## quartz_off_screen 
##                 2

Split

# split into 3 groups
hca.groups <- cutree(hca.ward, k=3)

# list of groups
print(sort(hca.groups))
##  CTC1 CTC10 CTC11 CTC12 CTC13 CTC15  CTC2  CTC3  CTC4  CTC5  CTC6  CTC7  CTC9 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##  NP29  NP35  NP36  NP37  NP38  NP39  NP41  NP42  NP43  NP44   S29   S37   S77 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##   S89 CTC14  NP20  NP30  NP34   S19   S21   S22   S27   S30   S31   S32   S33 
##     1     2     2     2     2     2     2     2     2     2     2     2     2 
##   S34   S35   S36   S39   S40   S41   S46   S47   S65  NP14   NP2  NP27   NP5 
##     2     2     2     2     2     2     2     2     2     3     3     3     3 
##   NP8  S100  S101  S102  S104  S105  S106  S107  S108  S109  S110  S113  S121 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##  S122  S123  S124  S126  S127  S128  S146  S147  S148  S150  S151  S152  S153 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##  S154  S160  S162  S163  S164  S165  S166  S167  S169  S170   S18   S20   S23 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S24   S25   S26   S28   S38   S42   S43   S44   S45   S48   S49   S50   S51 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S52   S53   S55   S56   S57   S58   S59   S60   S61   S63   S64   S79   S80 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S83   S84   S85   S86   S87   S88 
##     3     3     3     3     3     3

Make a dataframe from these groups

# df
wolbachia.groups <- print(sort(hca.groups)) %>% as.data.frame()
##  CTC1 CTC10 CTC11 CTC12 CTC13 CTC15  CTC2  CTC3  CTC4  CTC5  CTC6  CTC7  CTC9 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##  NP29  NP35  NP36  NP37  NP38  NP39  NP41  NP42  NP43  NP44   S29   S37   S77 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##   S89 CTC14  NP20  NP30  NP34   S19   S21   S22   S27   S30   S31   S32   S33 
##     1     2     2     2     2     2     2     2     2     2     2     2     2 
##   S34   S35   S36   S39   S40   S41   S46   S47   S65  NP14   NP2  NP27   NP5 
##     2     2     2     2     2     2     2     2     2     3     3     3     3 
##   NP8  S100  S101  S102  S104  S105  S106  S107  S108  S109  S110  S113  S121 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##  S122  S123  S124  S126  S127  S128  S146  S147  S148  S150  S151  S152  S153 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##  S154  S160  S162  S163  S164  S165  S166  S167  S169  S170   S18   S20   S23 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S24   S25   S26   S28   S38   S42   S43   S44   S45   S48   S49   S50   S51 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S52   S53   S55   S56   S57   S58   S59   S60   S61   S63   S64   S79   S80 
##     3     3     3     3     3     3     3     3     3     3     3     3     3 
##   S83   S84   S85   S86   S87   S88 
##     3     3     3     3     3     3
# change colnames
colnames(wolbachia.groups) <- "Group"

# change rownames
wolbachia.groups$Sample <- rownames(wolbachia.groups)

# add Group to otu_wolbachia with merge
wolbachia.groups <- wolbachia.groups %>% merge(otu_wolbachia %>% select(c(Sample, Wolbachia)), by="Sample")

# change rownames
rownames(wolbachia.groups) <- wolbachia.groups$Sample

# print
wolbachia.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Group Wolbachia
CTC1 CTC1 1 0.0000000
CTC10 CTC10 1 0.0352626
CTC11 CTC11 1 0.0077123
CTC12 CTC12 1 0.0191972
CTC13 CTC13 1 0.0026615
CTC14 CTC14 2 0.1074139
CTC15 CTC15 1 0.0094047
CTC2 CTC2 1 0.0087316
CTC3 CTC3 1 0.0007014
CTC4 CTC4 1 0.0088827
CTC5 CTC5 1 0.0017734
CTC6 CTC6 1 0.0020249
CTC7 CTC7 1 0.0050378
CTC9 CTC9 1 0.0058233
NP14 NP14 3 0.9888611
NP2 NP2 3 0.9993923
NP20 NP20 2 0.4705882
NP27 NP27 3 0.9789303
NP29 NP29 1 0.0492611
NP30 NP30 2 0.2894737
NP34 NP34 2 0.3157895
NP35 NP35 1 0.0032998
NP36 NP36 1 0.0000000
NP37 NP37 1 0.0002480
NP38 NP38 1 0.0218636
NP39 NP39 1 0.0023176
NP41 NP41 1 0.0000162
NP42 NP42 1 0.0004591
NP43 NP43 1 0.0070779
NP44 NP44 1 0.0000127
NP5 NP5 3 0.9999973
NP8 NP8 3 0.9999492
S100 S100 3 1.0000000
S101 S101 3 0.8433759
S102 S102 3 1.0000000
S104 S104 3 0.9999237
S105 S105 3 0.9999100
S106 S106 3 1.0000000
S107 S107 3 0.9999808
S108 S108 3 0.9999641
S109 S109 3 0.9999492
S110 S110 3 1.0000000
S113 S113 3 0.9941796
S121 S121 3 1.0000000
S122 S122 3 1.0000000
S123 S123 3 1.0000000
S124 S124 3 0.9999449
S126 S126 3 0.9997374
S127 S127 3 1.0000000
S128 S128 3 1.0000000
S146 S146 3 0.9993204
S147 S147 3 0.9999205
S148 S148 3 0.9993646
S150 S150 3 1.0000000
S151 S151 3 1.0000000
S152 S152 3 0.9998674
S153 S153 3 0.9998826
S154 S154 3 0.9993767
S160 S160 3 0.9999035
S162 S162 3 0.9999603
S163 S163 3 0.9999188
S164 S164 3 0.9993741
S165 S165 3 0.9995488
S166 S166 3 0.9992847
S167 S167 3 0.9997153
S169 S169 3 1.0000000
S170 S170 3 1.0000000
S18 S18 3 0.9909091
S19 S19 2 0.6928881
S20 S20 3 0.9376633
S21 S21 2 0.6642608
S22 S22 2 0.5062214
S23 S23 3 0.9421566
S24 S24 3 0.9054080
S25 S25 3 0.9970647
S26 S26 3 0.9733037
S27 S27 2 0.2155884
S28 S28 3 0.9646902
S29 S29 1 0.0098804
S30 S30 2 0.1662131
S31 S31 2 0.1847797
S32 S32 2 0.2115416
S33 S33 2 0.3796025
S34 S34 2 0.5882353
S35 S35 2 0.4507416
S36 S36 2 0.1426408
S37 S37 1 0.0420938
S38 S38 3 0.8371053
S39 S39 2 0.3274431
S40 S40 2 0.4336303
S41 S41 2 0.3043478
S42 S42 3 0.9997565
S43 S43 3 0.9369888
S44 S44 3 0.9795664
S45 S45 3 0.9815447
S46 S46 2 0.4921260
S47 S47 2 0.4081812
S48 S48 3 0.9938841
S49 S49 3 0.9925667
S50 S50 3 0.9898880
S51 S51 3 0.9936881
S52 S52 3 0.9833898
S53 S53 3 1.0000000
S55 S55 3 0.9999604
S56 S56 3 1.0000000
S57 S57 3 0.9998373
S58 S58 3 0.9977978
S59 S59 3 1.0000000
S60 S60 3 1.0000000
S61 S61 3 1.0000000
S63 S63 3 0.9955108
S64 S64 3 1.0000000
S65 S65 2 0.3235294
S77 S77 1 0.0000000
S79 S79 3 1.0000000
S80 S80 3 1.0000000
S83 S83 3 1.0000000
S84 S84 3 1.0000000
S85 S85 3 0.9994486
S86 S86 3 1.0000000
S87 S87 3 0.9998787
S88 S88 3 0.9977787
S89 S89 1 0.0032180

Groups statistics

With raw labels

# convert summary to table
summary.table <- function(df, n, var){
var <- df[df$Group==n, var]
x <- summary(var)
a <- data.frame(x=matrix(x),row.names=names(x))
colnames(a) <- n
return(a)
}

# make table with stats
wolbachia.summary <- summary.table(wolbachia.groups, 1, "Wolbachia") %>% cbind(summary.table(wolbachia.groups, 2, "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, 3, "Wolbachia"))

# print
wolbachia.summary %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
1 2 3
Min. 0.0000000 0.1074139 0.8371053
1st Qu. 0.0005802 0.2155884 0.9948452
Median 0.0032998 0.3274431 0.9998826
Mean 0.0091467 0.3654875 0.9891909
3rd Qu. 0.0091437 0.4705882 1.0000000
Max. 0.0492611 0.6928881 1.0000000

With new labels

# print levels of Group
levels(wolbachia.groups$Group %>% as.factor)
## [1] "1" "2" "3"
# change levels of Group
wolbachia.groups[wolbachia.groups$Group==1, "Group"] <- "Low infection"
wolbachia.groups[wolbachia.groups$Group==2, "Group"] <- "Medium infection"
wolbachia.groups[wolbachia.groups$Group==3, "Group"] <- "High infection"

# make table with stats
wolbachia.summary.labels <- summary.table(wolbachia.groups, "Low infection", "Wolbachia") %>% cbind(summary.table(wolbachia.groups, "Medium infection", "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, "High infection", "Wolbachia"))

# print
wolbachia.summary.labels %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Low infection Medium infection High infection
Min. 0.0000000 0.1074139 0.8371053
1st Qu. 0.0005802 0.2155884 0.9948452
Median 0.0032998 0.3274431 0.9998826
Mean 0.0091467 0.3654875 0.9891909
3rd Qu. 0.0091437 0.4705882 1.0000000
Max. 0.0492611 0.6928881 1.0000000

HCA on qPCR data

Preparation

# load qPCR file
setwd(path_qpcr)
qpcr <- read.xlsx("quanttiWolbachiaschriekeetal.xlsx")

# select Target/ref column
qpcr <- qpcr[,c("Sample.Name", "Target/Ref")]
colnames(qpcr) <- c("Sample", "Ratio_Target_Ref")

# remove the bad S45 sample
qpcr <- qpcr[-28,]
rownames(qpcr) <- qpcr$Sample

# samples without Wolbachia
qpcr_no_wolbachia <- qpcr[qpcr$Ratio_Target_Ref==0,]
qpcr_no_wolbachia$Group <- 0
qpcr_no_wolbachia <- qpcr_no_wolbachia %>% select(-c(Ratio_Target_Ref))

# samples with Wolbachia
qpcr_wolbachia <- qpcr[qpcr$Ratio_Target_Ref!=0,]

# print samples without Wolbachia
qpcr_no_wolbachia %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Group
CTC1 CTC1 0
CTC10 CTC10 0
CTC11 CTC11 0
CTC13 CTC13 0
CTC14 CTC14 0
CTC15 CTC15 0
CTC2 CTC2 0
CTC3 CTC3 0
CTC4 CTC4 0
CTC6 CTC6 0
CTC7 CTC7 0
# print samples with Wolbachia
qpcr_wolbachia %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Ratio_Target_Ref
CTC12 CTC12 0.000197
NJ39 NJ39 1.395000
NJ41 NJ41 3.481000
NJ42 NJ42 1.912000
NJ46 NJ46 0.496800
S18 S18 5.689000
S19 S19 13.310000
S20 S20 9.815000
S21 S21 4.573000
S22 S22 2.879000
S23 S23 7.873000
S25 S25 4.360000
S27 S27 0.114000
S30 S30 5.614000
S42 S42 2.207000
S43 S43 0.017600
S45 S45 0.352100
S47 S47 0.312800
S48 S48 0.237800
S50 S50 0.209600
S51 S51 3.615000
S52 S52 1.272000
S64 S64 2.427000
# add column with true ID in qpcr table
setwd(path_qpcr)
metadata <- openxlsx::read.xlsx("Supplementary_Table_1_V4_Shorten_w_ReadDepth_CheckqPCR.xlsx", sheet=1)
qpcr$True_TubeID <- rownames(qpcr)
qpcr <- qpcr %>% merge(metadata, by="True_TubeID", all.x=TRUE)
qpcr <- qpcr %>% select(c(True_TubeID, Sample.ID, Ratio_Target_Ref))

# print the qpcr table
qpcr %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
True_TubeID Sample.ID Ratio_Target_Ref
CTC1 CTC1 0.000000
CTC10 CTC10 0.000000
CTC11 CTC11 0.000000
CTC12 CTC12 0.000197
CTC13 CTC13 0.000000
CTC14 CTC14 0.000000
CTC15 CTC15 0.000000
CTC2 CTC2 0.000000
CTC3 CTC3 0.000000
CTC4 CTC4 0.000000
CTC6 CTC6 0.000000
CTC7 CTC7 0.000000
NJ39 NP27 1.395000
NJ41 NP29 3.481000
NJ42 NP30 1.912000
NJ46 NP34 0.496800
S18 S18 5.689000
S19 S19 13.310000
S20 S20 9.815000
S21 S21 4.573000
S22 S22 2.879000
S23 S23 7.873000
S25 S25 4.360000
S27 S27 0.114000
S30 S30 5.614000
S42 S42 2.207000
S43 S43 0.017600
S45 S45 0.352100
S47 S47 0.312800
S48 S48 0.237800
S50 S50 0.209600
S51 S51 3.615000
S52 S52 1.272000
S64 S64 2.427000

HCA

Run HCA

# Compute distance on samples with Wolbachia
d.wolbachia <- dist(qpcr_wolbachia[,-1, drop=FALSE], method="euclidean")

# HCA
hca.ward <- hclust(d.wolbachia,method="ward.D2")

Dendogram

# 3 groups
plot(hca.ward)
rect.hclust(hca.ward, k=3)

# 2 groups
plot(hca.ward)
rect.hclust(hca.ward,k=2)

Split

# split into 3 groups
hca.groups <- cutree(hca.ward, k=2)

# list of groups
print(sort(hca.groups))
## CTC12  NJ39  NJ41  NJ42  NJ46   S18   S21   S22   S25   S27   S30   S42   S43 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##   S45   S47   S48   S50   S51   S52   S64   S19   S20   S23 
##     1     1     1     1     1     1     1     2     2     2

Make a dataframe from these groups

Create dataframe

qpcr.groups <- print(sort(hca.groups)) %>% as.data.frame()
## CTC12  NJ39  NJ41  NJ42  NJ46   S18   S21   S22   S25   S27   S30   S42   S43 
##     1     1     1     1     1     1     1     1     1     1     1     1     1 
##   S45   S47   S48   S50   S51   S52   S64   S19   S20   S23 
##     1     1     1     1     1     1     1     2     2     2
colnames(qpcr.groups) <- "Group"
qpcr.groups$Sample <- rownames(qpcr.groups)

Add samples with no Wolbachia

qpcr.groups <- qpcr.groups %>% rbind(qpcr_no_wolbachia)

# print
qpcr.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Group Sample
CTC12 1 CTC12
NJ39 1 NJ39
NJ41 1 NJ41
NJ42 1 NJ42
NJ46 1 NJ46
S18 1 S18
S21 1 S21
S22 1 S22
S25 1 S25
S27 1 S27
S30 1 S30
S42 1 S42
S43 1 S43
S45 1 S45
S47 1 S47
S48 1 S48
S50 1 S50
S51 1 S51
S52 1 S52
S64 1 S64
S19 2 S19
S20 2 S20
S23 2 S23
CTC1 0 CTC1
CTC10 0 CTC10
CTC11 0 CTC11
CTC13 0 CTC13
CTC14 0 CTC14
CTC15 0 CTC15
CTC2 0 CTC2
CTC3 0 CTC3
CTC4 0 CTC4
CTC6 0 CTC6
CTC7 0 CTC7

Add groups in the qpcr table

rownames(qpcr) <- qpcr$True_TubeID
qpcr <- qpcr %>% merge(qpcr.groups, by="row.names")
rownames(qpcr) <- qpcr$Sample.ID
qpcr <- qpcr[,-1]
qpcr <- qpcr %>% select(-c(Sample))

# print the qpcr table
qpcr %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
True_TubeID Sample.ID Ratio_Target_Ref Group
CTC1 CTC1 CTC1 0.000000 0
CTC10 CTC10 CTC10 0.000000 0
CTC11 CTC11 CTC11 0.000000 0
CTC12 CTC12 CTC12 0.000197 1
CTC13 CTC13 CTC13 0.000000 0
CTC14 CTC14 CTC14 0.000000 0
CTC15 CTC15 CTC15 0.000000 0
CTC2 CTC2 CTC2 0.000000 0
CTC3 CTC3 CTC3 0.000000 0
CTC4 CTC4 CTC4 0.000000 0
CTC6 CTC6 CTC6 0.000000 0
CTC7 CTC7 CTC7 0.000000 0
NP27 NJ39 NP27 1.395000 1
NP29 NJ41 NP29 3.481000 1
NP30 NJ42 NP30 1.912000 1
NP34 NJ46 NP34 0.496800 1
S18 S18 S18 5.689000 1
S19 S19 S19 13.310000 2
S20 S20 S20 9.815000 2
S21 S21 S21 4.573000 1
S22 S22 S22 2.879000 1
S23 S23 S23 7.873000 2
S25 S25 S25 4.360000 1
S27 S27 S27 0.114000 1
S30 S30 S30 5.614000 1
S42 S42 S42 2.207000 1
S43 S43 S43 0.017600 1
S45 S45 S45 0.352100 1
S47 S47 S47 0.312800 1
S48 S48 S48 0.237800 1
S50 S50 S50 0.209600 1
S51 S51 S51 3.615000 1
S52 S52 S52 1.272000 1
S64 S64 S64 2.427000 1

Groups statistics

With raw labels

# make table with stats
qpcr.summary <- summary.table(qpcr, 0, "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, 1, "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, 2, "Ratio_Target_Ref"))

# print
qpcr.summary %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
0 1 2
Min. 0 0.000197 7.87300
1st Qu. 0 0.294050 8.84400
Median 0 1.653500 9.81500
Mean 0 2.058245 10.33267
3rd Qu. 0 3.514500 11.56250
Max. 0 5.689000 13.31000

With new labels

# print levels of Group
levels(qpcr$Group %>% as.factor)
## [1] "0" "1" "2"
# change levels of Group
qpcr[qpcr$Group==0, "Group"] <- "No infection"
qpcr[qpcr$Group==1, "Group"] <- "Low infection"
qpcr[qpcr$Group==2, "Group"] <- "High infection"

# make table with stats
qpcr.summary.labels <- summary.table(qpcr, "No infection", "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, "Low infection", "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, "High infection", "Ratio_Target_Ref"))

# print
qpcr.summary.labels %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
No infection Low infection High infection
Min. 0 0.000197 7.87300
1st Qu. 0 0.294050 8.84400
Median 0 1.653500 9.81500
Mean 0 2.058245 10.33267
3rd Qu. 0 3.514500 11.56250
Max. 0 5.689000 13.31000

NMDS (Non-metric MultiDimensional Scaling)

Preparation

Merge groups of data from phyloseq and data from qpcr

# print wolbachia.groups (phyloseq)
wolbachia.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Group Wolbachia
CTC1 CTC1 Low infection 0.0000000
CTC10 CTC10 Low infection 0.0352626
CTC11 CTC11 Low infection 0.0077123
CTC12 CTC12 Low infection 0.0191972
CTC13 CTC13 Low infection 0.0026615
CTC14 CTC14 Medium infection 0.1074139
CTC15 CTC15 Low infection 0.0094047
CTC2 CTC2 Low infection 0.0087316
CTC3 CTC3 Low infection 0.0007014
CTC4 CTC4 Low infection 0.0088827
CTC5 CTC5 Low infection 0.0017734
CTC6 CTC6 Low infection 0.0020249
CTC7 CTC7 Low infection 0.0050378
CTC9 CTC9 Low infection 0.0058233
NP14 NP14 High infection 0.9888611
NP2 NP2 High infection 0.9993923
NP20 NP20 Medium infection 0.4705882
NP27 NP27 High infection 0.9789303
NP29 NP29 Low infection 0.0492611
NP30 NP30 Medium infection 0.2894737
NP34 NP34 Medium infection 0.3157895
NP35 NP35 Low infection 0.0032998
NP36 NP36 Low infection 0.0000000
NP37 NP37 Low infection 0.0002480
NP38 NP38 Low infection 0.0218636
NP39 NP39 Low infection 0.0023176
NP41 NP41 Low infection 0.0000162
NP42 NP42 Low infection 0.0004591
NP43 NP43 Low infection 0.0070779
NP44 NP44 Low infection 0.0000127
NP5 NP5 High infection 0.9999973
NP8 NP8 High infection 0.9999492
S100 S100 High infection 1.0000000
S101 S101 High infection 0.8433759
S102 S102 High infection 1.0000000
S104 S104 High infection 0.9999237
S105 S105 High infection 0.9999100
S106 S106 High infection 1.0000000
S107 S107 High infection 0.9999808
S108 S108 High infection 0.9999641
S109 S109 High infection 0.9999492
S110 S110 High infection 1.0000000
S113 S113 High infection 0.9941796
S121 S121 High infection 1.0000000
S122 S122 High infection 1.0000000
S123 S123 High infection 1.0000000
S124 S124 High infection 0.9999449
S126 S126 High infection 0.9997374
S127 S127 High infection 1.0000000
S128 S128 High infection 1.0000000
S146 S146 High infection 0.9993204
S147 S147 High infection 0.9999205
S148 S148 High infection 0.9993646
S150 S150 High infection 1.0000000
S151 S151 High infection 1.0000000
S152 S152 High infection 0.9998674
S153 S153 High infection 0.9998826
S154 S154 High infection 0.9993767
S160 S160 High infection 0.9999035
S162 S162 High infection 0.9999603
S163 S163 High infection 0.9999188
S164 S164 High infection 0.9993741
S165 S165 High infection 0.9995488
S166 S166 High infection 0.9992847
S167 S167 High infection 0.9997153
S169 S169 High infection 1.0000000
S170 S170 High infection 1.0000000
S18 S18 High infection 0.9909091
S19 S19 Medium infection 0.6928881
S20 S20 High infection 0.9376633
S21 S21 Medium infection 0.6642608
S22 S22 Medium infection 0.5062214
S23 S23 High infection 0.9421566
S24 S24 High infection 0.9054080
S25 S25 High infection 0.9970647
S26 S26 High infection 0.9733037
S27 S27 Medium infection 0.2155884
S28 S28 High infection 0.9646902
S29 S29 Low infection 0.0098804
S30 S30 Medium infection 0.1662131
S31 S31 Medium infection 0.1847797
S32 S32 Medium infection 0.2115416
S33 S33 Medium infection 0.3796025
S34 S34 Medium infection 0.5882353
S35 S35 Medium infection 0.4507416
S36 S36 Medium infection 0.1426408
S37 S37 Low infection 0.0420938
S38 S38 High infection 0.8371053
S39 S39 Medium infection 0.3274431
S40 S40 Medium infection 0.4336303
S41 S41 Medium infection 0.3043478
S42 S42 High infection 0.9997565
S43 S43 High infection 0.9369888
S44 S44 High infection 0.9795664
S45 S45 High infection 0.9815447
S46 S46 Medium infection 0.4921260
S47 S47 Medium infection 0.4081812
S48 S48 High infection 0.9938841
S49 S49 High infection 0.9925667
S50 S50 High infection 0.9898880
S51 S51 High infection 0.9936881
S52 S52 High infection 0.9833898
S53 S53 High infection 1.0000000
S55 S55 High infection 0.9999604
S56 S56 High infection 1.0000000
S57 S57 High infection 0.9998373
S58 S58 High infection 0.9977978
S59 S59 High infection 1.0000000
S60 S60 High infection 1.0000000
S61 S61 High infection 1.0000000
S63 S63 High infection 0.9955108
S64 S64 High infection 1.0000000
S65 S65 Medium infection 0.3235294
S77 S77 Low infection 0.0000000
S79 S79 High infection 1.0000000
S80 S80 High infection 1.0000000
S83 S83 High infection 1.0000000
S84 S84 High infection 1.0000000
S85 S85 High infection 0.9994486
S86 S86 High infection 1.0000000
S87 S87 High infection 0.9998787
S88 S88 High infection 0.9977787
S89 S89 Low infection 0.0032180
# print qpcr.groups (qpcr)
qpcr %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
True_TubeID Sample.ID Ratio_Target_Ref Group
CTC1 CTC1 CTC1 0.000000 No infection
CTC10 CTC10 CTC10 0.000000 No infection
CTC11 CTC11 CTC11 0.000000 No infection
CTC12 CTC12 CTC12 0.000197 Low infection
CTC13 CTC13 CTC13 0.000000 No infection
CTC14 CTC14 CTC14 0.000000 No infection
CTC15 CTC15 CTC15 0.000000 No infection
CTC2 CTC2 CTC2 0.000000 No infection
CTC3 CTC3 CTC3 0.000000 No infection
CTC4 CTC4 CTC4 0.000000 No infection
CTC6 CTC6 CTC6 0.000000 No infection
CTC7 CTC7 CTC7 0.000000 No infection
NP27 NJ39 NP27 1.395000 Low infection
NP29 NJ41 NP29 3.481000 Low infection
NP30 NJ42 NP30 1.912000 Low infection
NP34 NJ46 NP34 0.496800 Low infection
S18 S18 S18 5.689000 Low infection
S19 S19 S19 13.310000 High infection
S20 S20 S20 9.815000 High infection
S21 S21 S21 4.573000 Low infection
S22 S22 S22 2.879000 Low infection
S23 S23 S23 7.873000 High infection
S25 S25 S25 4.360000 Low infection
S27 S27 S27 0.114000 Low infection
S30 S30 S30 5.614000 Low infection
S42 S42 S42 2.207000 Low infection
S43 S43 S43 0.017600 Low infection
S45 S45 S45 0.352100 Low infection
S47 S47 S47 0.312800 Low infection
S48 S48 S48 0.237800 Low infection
S50 S50 S50 0.209600 Low infection
S51 S51 S51 3.615000 Low infection
S52 S52 S52 1.272000 Low infection
S64 S64 S64 2.427000 Low infection
# merge 
colnames(wolbachia.groups) <- c("Sample", "Infection", "Wolbachia")
colnames(qpcr) <- c("qPCR_ID", "Sample", "qPCR_Ratio_Target_Ref", "qPCR_Infection")
all.groups <- wolbachia.groups %>% merge(qpcr, by="Sample", all.x=TRUE)

all.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample Infection Wolbachia qPCR_ID qPCR_Ratio_Target_Ref qPCR_Infection
CTC1 Low infection 0.0000000 CTC1 0.000000 No infection
CTC10 Low infection 0.0352626 CTC10 0.000000 No infection
CTC11 Low infection 0.0077123 CTC11 0.000000 No infection
CTC12 Low infection 0.0191972 CTC12 0.000197 Low infection
CTC13 Low infection 0.0026615 CTC13 0.000000 No infection
CTC14 Medium infection 0.1074139 CTC14 0.000000 No infection
CTC15 Low infection 0.0094047 CTC15 0.000000 No infection
CTC2 Low infection 0.0087316 CTC2 0.000000 No infection
CTC3 Low infection 0.0007014 CTC3 0.000000 No infection
CTC4 Low infection 0.0088827 CTC4 0.000000 No infection
CTC5 Low infection 0.0017734 NA NA NA
CTC6 Low infection 0.0020249 CTC6 0.000000 No infection
CTC7 Low infection 0.0050378 CTC7 0.000000 No infection
CTC9 Low infection 0.0058233 NA NA NA
NP14 High infection 0.9888611 NA NA NA
NP2 High infection 0.9993923 NA NA NA
NP20 Medium infection 0.4705882 NA NA NA
NP27 High infection 0.9789303 NJ39 1.395000 Low infection
NP29 Low infection 0.0492611 NJ41 3.481000 Low infection
NP30 Medium infection 0.2894737 NJ42 1.912000 Low infection
NP34 Medium infection 0.3157895 NJ46 0.496800 Low infection
NP35 Low infection 0.0032998 NA NA NA
NP36 Low infection 0.0000000 NA NA NA
NP37 Low infection 0.0002480 NA NA NA
NP38 Low infection 0.0218636 NA NA NA
NP39 Low infection 0.0023176 NA NA NA
NP41 Low infection 0.0000162 NA NA NA
NP42 Low infection 0.0004591 NA NA NA
NP43 Low infection 0.0070779 NA NA NA
NP44 Low infection 0.0000127 NA NA NA
NP5 High infection 0.9999973 NA NA NA
NP8 High infection 0.9999492 NA NA NA
S100 High infection 1.0000000 NA NA NA
S101 High infection 0.8433759 NA NA NA
S102 High infection 1.0000000 NA NA NA
S104 High infection 0.9999237 NA NA NA
S105 High infection 0.9999100 NA NA NA
S106 High infection 1.0000000 NA NA NA
S107 High infection 0.9999808 NA NA NA
S108 High infection 0.9999641 NA NA NA
S109 High infection 0.9999492 NA NA NA
S110 High infection 1.0000000 NA NA NA
S113 High infection 0.9941796 NA NA NA
S121 High infection 1.0000000 NA NA NA
S122 High infection 1.0000000 NA NA NA
S123 High infection 1.0000000 NA NA NA
S124 High infection 0.9999449 NA NA NA
S126 High infection 0.9997374 NA NA NA
S127 High infection 1.0000000 NA NA NA
S128 High infection 1.0000000 NA NA NA
S146 High infection 0.9993204 NA NA NA
S147 High infection 0.9999205 NA NA NA
S148 High infection 0.9993646 NA NA NA
S150 High infection 1.0000000 NA NA NA
S151 High infection 1.0000000 NA NA NA
S152 High infection 0.9998674 NA NA NA
S153 High infection 0.9998826 NA NA NA
S154 High infection 0.9993767 NA NA NA
S160 High infection 0.9999035 NA NA NA
S162 High infection 0.9999603 NA NA NA
S163 High infection 0.9999188 NA NA NA
S164 High infection 0.9993741 NA NA NA
S165 High infection 0.9995488 NA NA NA
S166 High infection 0.9992847 NA NA NA
S167 High infection 0.9997153 NA NA NA
S169 High infection 1.0000000 NA NA NA
S170 High infection 1.0000000 NA NA NA
S18 High infection 0.9909091 S18 5.689000 Low infection
S19 Medium infection 0.6928881 S19 13.310000 High infection
S20 High infection 0.9376633 S20 9.815000 High infection
S21 Medium infection 0.6642608 S21 4.573000 Low infection
S22 Medium infection 0.5062214 S22 2.879000 Low infection
S23 High infection 0.9421566 S23 7.873000 High infection
S24 High infection 0.9054080 NA NA NA
S25 High infection 0.9970647 S25 4.360000 Low infection
S26 High infection 0.9733037 NA NA NA
S27 Medium infection 0.2155884 S27 0.114000 Low infection
S28 High infection 0.9646902 NA NA NA
S29 Low infection 0.0098804 NA NA NA
S30 Medium infection 0.1662131 S30 5.614000 Low infection
S31 Medium infection 0.1847797 NA NA NA
S32 Medium infection 0.2115416 NA NA NA
S33 Medium infection 0.3796025 NA NA NA
S34 Medium infection 0.5882353 NA NA NA
S35 Medium infection 0.4507416 NA NA NA
S36 Medium infection 0.1426408 NA NA NA
S37 Low infection 0.0420938 NA NA NA
S38 High infection 0.8371053 NA NA NA
S39 Medium infection 0.3274431 NA NA NA
S40 Medium infection 0.4336303 NA NA NA
S41 Medium infection 0.3043478 NA NA NA
S42 High infection 0.9997565 S42 2.207000 Low infection
S43 High infection 0.9369888 S43 0.017600 Low infection
S44 High infection 0.9795664 NA NA NA
S45 High infection 0.9815447 S45 0.352100 Low infection
S46 Medium infection 0.4921260 NA NA NA
S47 Medium infection 0.4081812 S47 0.312800 Low infection
S48 High infection 0.9938841 S48 0.237800 Low infection
S49 High infection 0.9925667 NA NA NA
S50 High infection 0.9898880 S50 0.209600 Low infection
S51 High infection 0.9936881 S51 3.615000 Low infection
S52 High infection 0.9833898 S52 1.272000 Low infection
S53 High infection 1.0000000 NA NA NA
S55 High infection 0.9999604 NA NA NA
S56 High infection 1.0000000 NA NA NA
S57 High infection 0.9998373 NA NA NA
S58 High infection 0.9977978 NA NA NA
S59 High infection 1.0000000 NA NA NA
S60 High infection 1.0000000 NA NA NA
S61 High infection 1.0000000 NA NA NA
S63 High infection 0.9955108 NA NA NA
S64 High infection 1.0000000 S64 2.427000 Low infection
S65 Medium infection 0.3235294 NA NA NA
S77 Low infection 0.0000000 NA NA NA
S79 High infection 1.0000000 NA NA NA
S80 High infection 1.0000000 NA NA NA
S83 High infection 1.0000000 NA NA NA
S84 High infection 1.0000000 NA NA NA
S85 High infection 0.9994486 NA NA NA
S86 High infection 1.0000000 NA NA NA
S87 High infection 0.9998787 NA NA NA
S88 High infection 0.9977787 NA NA NA
S89 Low infection 0.0032180 NA NA NA
# optionnal : all.groups for comparison
all.groups.comparison <- wolbachia.groups %>% merge(qpcr, by="Sample")

all.groups.comparison <- all.groups.comparison %>% select(c(Sample, qPCR_ID, Infection, Wolbachia, qPCR_Ratio_Target_Ref, qPCR_Infection))

all.groups.comparison$Infection <- all.groups.comparison$Infection %>% as.factor()

all.groups.comparison$Infection <- factor(all.groups.comparison$Infection, levels = c("Low infection", "Medium infection", "High infection"))

all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]

all.groups.comparison %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample qPCR_ID Infection Wolbachia qPCR_Ratio_Target_Ref qPCR_Infection
1 CTC1 CTC1 Low infection 0.0000000 0.000000 No infection
2 CTC10 CTC10 Low infection 0.0352626 0.000000 No infection
3 CTC11 CTC11 Low infection 0.0077123 0.000000 No infection
4 CTC12 CTC12 Low infection 0.0191972 0.000197 Low infection
5 CTC13 CTC13 Low infection 0.0026615 0.000000 No infection
7 CTC15 CTC15 Low infection 0.0094047 0.000000 No infection
8 CTC2 CTC2 Low infection 0.0087316 0.000000 No infection
9 CTC3 CTC3 Low infection 0.0007014 0.000000 No infection
10 CTC4 CTC4 Low infection 0.0088827 0.000000 No infection
11 CTC6 CTC6 Low infection 0.0020249 0.000000 No infection
12 CTC7 CTC7 Low infection 0.0050378 0.000000 No infection
14 NP29 NJ41 Low infection 0.0492611 3.481000 Low infection
6 CTC14 CTC14 Medium infection 0.1074139 0.000000 No infection
15 NP30 NJ42 Medium infection 0.2894737 1.912000 Low infection
16 NP34 NJ46 Medium infection 0.3157895 0.496800 Low infection
18 S19 S19 Medium infection 0.6928881 13.310000 High infection
20 S21 S21 Medium infection 0.6642608 4.573000 Low infection
21 S22 S22 Medium infection 0.5062214 2.879000 Low infection
24 S27 S27 Medium infection 0.2155884 0.114000 Low infection
25 S30 S30 Medium infection 0.1662131 5.614000 Low infection
29 S47 S47 Medium infection 0.4081812 0.312800 Low infection
13 NP27 NJ39 High infection 0.9789303 1.395000 Low infection
17 S18 S18 High infection 0.9909091 5.689000 Low infection
19 S20 S20 High infection 0.9376633 9.815000 High infection
22 S23 S23 High infection 0.9421566 7.873000 High infection
23 S25 S25 High infection 0.9970647 4.360000 Low infection
26 S42 S42 High infection 0.9997565 2.207000 Low infection
27 S43 S43 High infection 0.9369888 0.017600 Low infection
28 S45 S45 High infection 0.9815447 0.352100 Low infection
30 S48 S48 High infection 0.9938841 0.237800 Low infection
31 S50 S50 High infection 0.9898880 0.209600 Low infection
32 S51 S51 High infection 0.9936881 3.615000 Low infection
33 S52 S52 High infection 0.9833898 1.272000 Low infection
34 S64 S64 High infection 1.0000000 2.427000 Low infection
all.groups.comparison <- all.groups.comparison %>% merge(sam %>% select(c(Sample, Strain, Species)), by="Sample")

all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]

all.groups.comparison %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)
Sample qPCR_ID Infection Wolbachia qPCR_Ratio_Target_Ref qPCR_Infection Strain Species
1 CTC1 CTC1 Low infection 0.0000000 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
2 CTC10 CTC10 Low infection 0.0352626 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
3 CTC11 CTC11 Low infection 0.0077123 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
4 CTC12 CTC12 Low infection 0.0191972 0.000197 Low infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
5 CTC13 CTC13 Low infection 0.0026615 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
7 CTC15 CTC15 Low infection 0.0094047 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
8 CTC2 CTC2 Low infection 0.0087316 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
9 CTC3 CTC3 Low infection 0.0007014 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
10 CTC4 CTC4 Low infection 0.0088827 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
11 CTC6 CTC6 Low infection 0.0020249 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
12 CTC7 CTC7 Low infection 0.0050378 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
14 NP29 NJ41 Low infection 0.0492611 3.481000 Low infection Field - Guadeloupe Culex quinquefasciatus
6 CTC14 CTC14 Medium infection 0.1074139 0.000000 No infection Laboratory - Slab TC (Wolbachia -) Culex quinquefasciatus
15 NP30 NJ42 Medium infection 0.2894737 1.912000 Low infection Field - Guadeloupe Culex quinquefasciatus
16 NP34 NJ46 Medium infection 0.3157895 0.496800 Low infection Field - Guadeloupe Culex quinquefasciatus
18 S19 S19 Medium infection 0.6928881 13.310000 High infection Laboratory - Lavar Culex pipiens
20 S21 S21 Medium infection 0.6642608 4.573000 Low infection Laboratory - Lavar Culex pipiens
21 S22 S22 Medium infection 0.5062214 2.879000 Low infection Laboratory - Lavar Culex pipiens
24 S27 S27 Medium infection 0.2155884 0.114000 Low infection Laboratory - Lavar Culex pipiens
25 S30 S30 Medium infection 0.1662131 5.614000 Low infection Laboratory - Lavar Culex pipiens
29 S47 S47 Medium infection 0.4081812 0.312800 Low infection Field - Camping Europe Culex pipiens
13 NP27 NJ39 High infection 0.9789303 1.395000 Low infection Field - Guadeloupe Culex quinquefasciatus
17 S18 S18 High infection 0.9909091 5.689000 Low infection Laboratory - Lavar Culex pipiens
19 S20 S20 High infection 0.9376633 9.815000 High infection Laboratory - Lavar Culex pipiens
22 S23 S23 High infection 0.9421566 7.873000 High infection Laboratory - Lavar Culex pipiens
23 S25 S25 High infection 0.9970647 4.360000 Low infection Laboratory - Lavar Culex pipiens
26 S42 S42 High infection 0.9997565 2.207000 Low infection Field - Camping Europe Culex pipiens
27 S43 S43 High infection 0.9369888 0.017600 Low infection Field - Camping Europe Culex pipiens
28 S45 S45 High infection 0.9815447 0.352100 Low infection Field - Camping Europe Culex pipiens
30 S48 S48 High infection 0.9938841 0.237800 Low infection Field - Camping Europe Culex pipiens
31 S50 S50 High infection 0.9898880 0.209600 Low infection Field - Bosc Culex pipiens
32 S51 S51 High infection 0.9936881 3.615000 Low infection Field - Bosc Culex pipiens
33 S52 S52 High infection 0.9833898 1.272000 Low infection Field - Bosc Culex pipiens
34 S64 S64 High infection 1.0000000 2.427000 Low infection Field - Bosc Culex pipiens
means <- aggregate(qPCR_Ratio_Target_Ref ~ Infection , all.groups.comparison, median)
means$qPCR_Ratio_Target_Ref <- means$qPCR_Ratio_Target_Ref %>% round(2)

ggqqplot(all.groups.comparison$Wolbachia)

shapiro.test(all.groups.comparison$Wolbachia)
## 
##  Shapiro-Wilk normality test
## 
## data:  all.groups.comparison$Wolbachia
## W = 0.78016, p-value = 1.084e-05
my_comparisons <- list( c("Low infection", "Medium infection"), 
                        c("Medium infection", "High infection"), 
                        c("Low infection", "High infection") )

p1 <- ggplot(all.groups.comparison, aes(x=Infection, y=qPCR_Ratio_Target_Ref, fill=Infection)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
        geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
        theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")), 
       y="Ratio Target Ref (qPCR)", 
       fill=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")))+
  guides(colour=FALSE)+
  #stat_summary(fun.y = "median", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(10, 11, 12))+
    stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(10, 11, 12))+
  geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.6))

p1

p9 <- ggplot(all.groups.comparison, aes(x=Infection, y=Wolbachia, fill=Infection)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")), 
       y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")), 
       fill="Infection")+
  guides(colour=FALSE)+ 
  scale_y_continuous(breaks = seq(0, 1, by = 0.1))

p9

setwd(path_tsv)
write.table(all.groups.comparison, "1Ec_qpcr_comparison.tsv", sep="\t", row.names = FALSE, col.names=TRUE)

setwd(path_plot)
tiff("1Ec_boxplot.tiff", units="in", width=10, height=5, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2
png("1Ec_boxplot.png", units="in", width=10, height=5, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2

Add groups in decontam phyloseq object

# load phyloseq object after decontam
setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")

# extract metadata
sam <- data.frame(ps.filter@sam_data)

# add groups in sam table
sam <- sam %>% merge(all.groups %>% select(c(Sample, Infection, Wolbachia, qPCR_Infection, qPCR_Ratio_Target_Ref)), by="Sample", all.x = TRUE)
rownames(sam)<- sam$Sample

# change levels for plot
# levels(sam$Strain) <- c(levels(sam$Strain), "Laboratory - Lavar", "Laboratory - Slab TC (Wolbachia -)")
# sam$Strain[sam$Strain=="Laboratory - Lavar"] <- "Laboratory - Lavar"
# sam$Strain[sam$Strain=="Wolbachia -"] <- "Laboratory - Slab TC (Wolbachia -)"
# sam$Strain <- droplevels(sam$Strain)
# levels(sam$Strain)

# add a column for qPCR measure
sam$qPCR_measure <- "Yes"
sam[is.na(sam$qPCR_Ratio_Target_Ref), "qPCR_measure"] <- "No"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref==0, "qPCR_measure"] <- "= 0 (Wolbachia -)"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref>0, "qPCR_measure"] <- "> 0 (Wolbachia +)"

# update sam table of phyloseq object
sample_data(ps.filter) <- sam

# transformation percent
ps_percent <- transform_sample_counts(ps.filter, function(x) x/sum(x)*100)

# select whole
ps_percent_whole <- subset_samples(ps_percent, Organ == "Whole")

Compute Bray-Curtis distance

set.seed(35)
bray.culex_whole <- ordinate(ps_percent_whole, method="NMDS", distance="bray", trymax=300)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1088607 
## Run 1 stress 0.1070017 
## ... New best solution
## ... Procrustes: rmse 0.02445929  max resid 0.1700959 
## Run 2 stress 0.1452511 
## Run 3 stress 0.1087863 
## Run 4 stress 0.1522622 
## Run 5 stress 0.1215459 
## Run 6 stress 0.1553544 
## Run 7 stress 0.119003 
## Run 8 stress 0.1259113 
## Run 9 stress 0.1064512 
## ... New best solution
## ... Procrustes: rmse 0.006708485  max resid 0.03812992 
## Run 10 stress 0.1397327 
## Run 11 stress 0.1215459 
## Run 12 stress 0.1064512 
## ... Procrustes: rmse 1.293125e-05  max resid 8.41373e-05 
## ... Similar to previous best
## Run 13 stress 0.1521948 
## Run 14 stress 0.1389989 
## Run 15 stress 0.1558452 
## Run 16 stress 0.1090684 
## Run 17 stress 0.1174935 
## Run 18 stress 0.1090611 
## Run 19 stress 0.1070004 
## Run 20 stress 0.1175567 
## *** Solution reached

Plot

Raw data

# all info
p1 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape=NA) +
  labs(x="NMDS1", y = "NMDS2")+
  stat_ellipse(geom = "polygon", level=0.95, alpha = 0.25, aes(fill = Species, color=NA))+
  scale_fill_manual(values=c("orange", "yellow","green"), labels = c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
  scale_color_manual(values=c("blue", "grey","red"), breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_shape_manual(values=c(19,17,15,4,7), labels=c("Field - Bosc", "Field - Camping Europe", "Field - Guadeloupe", "Laboratory - Lavar", expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)"))))+
  theme_gray()+
  labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  guides(fill = guide_legend(order = 1),
         colour = guide_legend(order = 2),
         shape = guide_legend(order = 3))+
  labs(fill="Species", col="Wolbachia Group")+
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)
## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "Infection", : Shape variable was not found in the available data you
## provided.No shape mapped.
p1$layers[[1]] <- NULL
p1[["guides"]][["shape"]][["title"]] <- "Strain"
p1[["guides"]][["colour"]][["title"]] <- expression(paste(italic("Wolbachia"), "Infection"))

p1+
  geom_point(aes(shape = Strain), size = 2)

# Wolbachia infection

p2 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape="qPCR_measure") +
  
  scale_color_manual(values=c("blue", "#E4A7F5","red"), 
                     breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_fill_manual(expression(paste(italic("Wolbachia"), "Infection")), 
                    values=c("blue", "#E4A7F5","red"), 
                    breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_shape_manual("qPCR", 
                     values=c(7,2,19), 
                     labels=c(expression(paste("= 0 (", italic("Wolbachia"), "-)")), 
                              expression(paste("> 0 (", italic("Wolbachia"), "+)")), 
                              "No qPCR"))+
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  
  stat_ellipse(geom = "polygon", 
               alpha = 0.25, 
               aes(fill=Infection, group=Infection))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         shape = guide_legend(order = 2),
         colour=FALSE)+
  
  theme(legend.text.align = 0)

p2

# species ("orange", "yellow", "green")
p3 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Species") +
  
  geom_point(size = 2) +
  
  scale_fill_manual(values=c("#F57462", "#498CF5","#559E36"), 
                    labels = c(expression(italic("Aedes aegypti")), 
                               expression(italic("Culex pipiens")), 
                               expression(italic("Culex quinquefasciatus"))))+
  
  scale_color_manual(values=c("#F57462", "#498CF5","#559E36"))+
  
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  
  stat_ellipse(geom = "polygon", 
               alpha = 0.25, 
               aes(fill = Species))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  
  theme(legend.text.align = 0)

p3

# Strain
new_location <- c("Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)", 
                  "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar")
ps_percent_whole@sam_data$Strain <- factor(ps_percent_whole@sam_data$Strain, levels = new_location)

p4 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain") +
  
  geom_point(size = 2) +
  
  scale_fill_manual("Strain", 
                    values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                    labels = c("Field - Guadeloupe", 
                               expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                               "Field - Bosc", 
                               "Field - Camping Europe", 
                               "Laboratory - Lavar"))+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                                "Field - Bosc", 
                                "Field - Camping Europe", 
                                "Laboratory - Lavar"))+
  
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Strain, group=Strain))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)

p4

# field
p5 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Field") +
  geom_point(size = 2) +
  scale_fill_manual(values=c("brown", "purple" ))+
  scale_color_manual(values=c("brown", "purple"))+
  
  labs(x="NMDS1", 
       y = "NMDS2", 
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Field, group=Field))+
  
  #coord_cartesian(xlim = c(-2,2.8), ylim=c(-2,1.5))+
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)

p5

# species + strains
p6 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain", shape=NA) +
  labs(x="NMDS1", y = "NMDS2")+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                                "Field - Bosc", 
                                "Field - Camping Europe", 
                                "Laboratory - Lavar"))+
  scale_shape_manual(values = c(15,16,17), labels=c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
  theme_gray()+
  labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  guides(fill = guide_legend(order = 1),
         colour = guide_legend(order = 2),
         shape = guide_legend(order = 3))+
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)+
  geom_point(aes(shape = Species_italic), size = 2)
## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "Strain", : Shape variable was not found in the available data you provided.No
## shape mapped.
p6$layers[[1]] <- NULL
p6[["guides"]][["shape"]][["title"]] <- "Species"

p6

sam2 <- sam[sam$Species=="Culex pipiens",]

means <- aggregate(Wolbachia ~ Strain , sam2, median)
means$Wolbachia <- means$Wolbachia %>% round(2)

ggqqplot(sam2$Wolbachia)

shapiro.test(sam2$Wolbachia)
## 
##  Shapiro-Wilk normality test
## 
## data:  sam2$Wolbachia
## W = 0.53983, p-value = 1.266e-14
my_comparisons <- list( c("Field - Bosc", "Field - Camping Europe"), 
                        c("Field - Camping Europe", "Laboratory - Lavar"), 
                        c("Field - Bosc", "Laboratory - Lavar") )

p7 <- ggplot(sam2, aes(x=Strain, y=Wolbachia, fill=Strain)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x="\nStrain", 
       y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")),
       fill="Strain")+
  guides(colour=FALSE)+
  #stat_summary(fun = "mean", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(1.1, 1.2, 1.3))+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(1.1, 1.2, 1.3))+
  geom_text(data = means, aes(label = Wolbachia, y = Wolbachia - 0.05))

p7
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.0000000002, 1.0000000002, 0.99992370072, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(0.9999999998, 1.0000000001, 1.0000000001, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(1.0000000002, 1.0000000002, 0.99992370072, :
## cannot compute exact p-value with ties

means <- aggregate(qPCR_Ratio_Target_Ref ~ Strain , sam2, median)
means$qPCR_Ratio_Target_Ref<- means$qPCR_Ratio_Target_Ref %>% round(2)

wilcox.test(qPCR_Ratio_Target_Ref ~ Strain, sam2[sam2$Strain!="Field - Camping Europe",])
## 
##  Wilcoxon rank sum test
## 
## data:  qPCR_Ratio_Target_Ref by Strain
## W = 5, p-value = 0.05035
## alternative hypothesis: true location shift is not equal to 0
sam3 <- na.omit(sam2)
ggqqplot(sam3$qPCR_Ratio_Target_Ref)

shapiro.test(sam3$qPCR_Ratio_Target_Ref)
## 
##  Shapiro-Wilk normality test
## 
## data:  sam3$qPCR_Ratio_Target_Ref
## W = 0.86499, p-value = 0.01473
p8 <- ggplot(na.omit(sam2), aes(x=Strain, y=qPCR_Ratio_Target_Ref, fill=Strain)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x="\nStrain", y="Ratio Target Ref (qPCR)", fill="Strain")+
  guides(colour=FALSE)+
  #stat_summary(fun = "median", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(10, 11, 12))+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(10, 11, 12))+
  geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.4))

p8

ggplot(na.omit(sam2), aes(x=Wolbachia, y=qPCR_Ratio_Target_Ref))+
  geom_line()

x <- sam3$Wolbachia %>% log10
y <- sam3$qPCR_Ratio_Target_Ref %>% log10

plot(x, y, pch = 19, col = "lightblue")
abline(lm(y ~ x), col = "red", lwd = 3)
text(paste("Correlation:", round(cor(x, y), 2)), x = -0.6, y = 0.2)

# panels
p_group <- plot_grid(p2,
                     p3,
                     p4,
                     p5,
                     nrow=2, 
                     ncol=2,
                     align="hv",
                     labels=c("A","B","C", "D")
                     )

p_group

p_group2 <- plot_grid(p7,
                      p8,
                      nrow=1, 
                      align="h",
                     labels=c("A","B"))

p_group2

#ggarrange(p2, p3, p4, p5,nrow=2)

# save
setwd(path_plot)

tiff("1Ec_NMDS_main.tiff", units="in", width=8, height=5, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ec_NMDS_supp.tiff", units="in", width=16, height=5, res=300)
p_group2
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ec_MED_panels.tiff", units="in", width=14, height=10, res=300)
p_group
dev.off()
## quartz_off_screen 
##                 2
png("1Ec_NMDS_main.png", units="in", width=8, height=5, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
png("1Ec_NMDS_supp.png", units="in", width=16, height=5, res=300)
p_group2
dev.off()
## quartz_off_screen 
##                 2
png("1Ec_MED_panels.png", units="in", width=14, height=10, res=300)
p_group
dev.off()
## quartz_off_screen 
##                 2

qPCR data

# all info
p1 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="qPCR_Infection", shape=NA) +
  labs(x="NMDS1", y = "NMDS2")+
  stat_ellipse(geom = "polygon", level=0.95, alpha = 0.25, aes(fill = Species, color=NA))+
  scale_fill_manual(values=c("orange", "yellow","green"), labels = c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
  scale_color_manual(values=c("black", "blue", "red"), breaks=c("No infection", "Low infection", "High infection"))+
  theme_gray()+
  labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  guides(fill = guide_legend(order = 1),
         colour = guide_legend(order = 2),
         shape = guide_legend(order = 3))+
  labs(fill="Species", col="Wolbachia Group")+
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)
## Warning in plot_ordination(ps_percent_whole, bray.culex_whole, color =
## "qPCR_Infection", : Shape variable was not found in the available data you
## provided.No shape mapped.
p1$layers[[1]] <- NULL
p1[["guides"]][["shape"]][["title"]] <- "Strain"
p1[["guides"]][["colour"]][["title"]] <- expression(paste(italic("Wolbachia"), "Infection"))

p1+
  geom_point(aes(shape = Strain), size = 2)
## Warning: Removed 34 rows containing missing values (geom_point).

Statistical analysis

Effect of field

require(vegan)
setwd(path_xlsx)

# Select pipiens whole from 
ps_pipiens <- ps.filter %>% subset_samples(Species=="Culex pipiens")
ps_pipiens <- check_ps(ps_pipiens)
ps_pipiens_whole <- ps_pipiens %>% subset_samples(Organ=="Whole")
ps_pipiens_whole_bosc_lavar <- subset_samples(ps_pipiens_whole, Strain!="Field - Camping Europe")
ps_pipiens_whole_camping_lavar <- subset_samples(ps_pipiens_whole, Strain!="Field - Bosc")
ps_pipiens_whole_bosc_camping <- subset_samples(ps_pipiens_whole, Strain!="Laboratory - Lavar")

# Adonis method

## Bosc / Lavar
adonis_field_bosc_lavar <- adonis(vegdist(t(otu_table(ps_pipiens_whole_bosc_lavar)), method = "bray") ~Field,
       data=as(sample_data(ps_pipiens_whole_bosc_lavar), "data.frame"), permutation = 9999)

res_adonis1 <- adonis_field_bosc_lavar$aov.tab
res_adonis1 <- fill_significance(res_adonis1,"Pr(>F)")

## Camping / Lavar
adonis_field_camping_lavar <- adonis(vegdist(t(otu_table(ps_pipiens_whole_camping_lavar)), method = "bray") ~Field,
       data=as(sample_data(ps_pipiens_whole_camping_lavar), "data.frame"), permutation = 9999)

res_adonis2 <- adonis_field_camping_lavar$aov.tab
res_adonis2 <- fill_significance(res_adonis2,"Pr(>F)")

## Bosc / Camping
adonis_field_bosc_camping <- adonis(vegdist(t(otu_table(ps_pipiens_whole_bosc_camping)), method = "bray") ~Strain,
       data=as(sample_data(ps_pipiens_whole_bosc_camping), "data.frame"), permutation = 9999)

res_adonis3 <- adonis_field_bosc_camping$aov.tab
res_adonis3 <- fill_significance(res_adonis3,"Pr(>F)")

# Excel
wb <- loadWorkbook("Supplementary_Table_2.xlsx")
addWorksheet(wb, "NMDS")

writeData(wb, "NMDS", "Effect of lab (and Strain) in the Culex pipiens samples", startCol = 1, startRow = 25)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 27)
writeData(wb, "NMDS", "Lavar / Bosc:", startCol = 1, startRow = 29)
writeData(wb, "NMDS", res_adonis1, startCol = 1, startRow = 30, rowNames=TRUE)
writeData(wb, "NMDS", "Lavar / Camping Europe:", startCol = 1, startRow = 35)
writeData(wb, "NMDS", res_adonis2, startCol = 1, startRow = 36, rowNames=TRUE)
writeData(wb, "NMDS", "Bosc / Camping Europe:", startCol = 1, startRow = 41)
writeData(wb, "NMDS", res_adonis3, startCol = 1, startRow = 42, rowNames=TRUE)


saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of tetracycline

setwd(path_xlsx)

# Select quinque whole
ps_quinque <- ps.filter %>% subset_samples(Species=="Culex quinquefasciatus")
ps_quinque <- check_ps(ps_quinque)
ps_quinque_whole <- ps_quinque %>% subset_samples(Organ=="Whole")

# Adonis method
adonis_antibiotics <- adonis(vegdist(t(otu_table(ps_quinque_whole)), method = "bray") ~Field,
       data=as(sample_data(ps_quinque_whole), "data.frame"), permutation = 9999)

adonis_antibiotics$aov.tab 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      1    1.4603 1.46029   5.191 0.25709  4e-04 ***
## Residuals 15    4.2197 0.28131         0.74291           
## Total     16    5.6800                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res_adonis5 <- adonis_antibiotics$aov.tab %>% as.data.frame()
res_adonis5 <- fill_significance(res_adonis3,"Pr(>F)")

writeData(wb, "NMDS", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 13)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 15)
writeData(wb, "NMDS", res_adonis3, startCol = 1, startRow = 17, rowNames=TRUE)

saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of species

setwd(path_xlsx)

ps_whole_field <- ps.filter %>% subset_samples(Field=="Field")

adonis_whole_field <- adonis(vegdist(t(otu_table(ps_whole_field)), method = "bray") ~Species,
       data=as(sample_data(ps_whole_field), "data.frame"), permutation = 9999)

res_adonis6 <- adonis_whole_field$aov.tab
res_adonis6 <- fill_significance(res_adonis6,"Pr(>F)")

writeData(wb, "NMDS", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 3)
writeData(wb, "NMDS", res_adonis6, startCol = 1, startRow = 5, rowNames=TRUE)

saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of Wolbachia infection within Culex pipiens

setwd(path_xlsx)

# HCA Wolbachia infection * Strain
adonis_whole_infection <- adonis(vegdist(t(otu_table(ps_pipiens_whole)), method = "bray") ~Infection*Strain,
       data=as(sample_data(ps_pipiens_whole), "data.frame"), permutation = 9999)

res_adonis7 <- adonis_whole_infection$aov.tab
res_adonis7 <- fill_significance(res_adonis7,"Pr(>F)")

writeData(wb, "NMDS", "Effect of Wolbachia infection (HCA groups) and Strain in Culex pipiens", startCol = 1, startRow = 50)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 52)
writeData(wb, "NMDS", res_adonis7, startCol = 1, startRow = 54, rowNames=TRUE)

saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of Wolbachia infection, qPCR values and Strain within Culex pipiens with qPCR values

setwd(path_xlsx)

# HCA Wolbachia infection * Strain in Culex pipiens with a qPCR value
ps_pipiens_whole_qpcr <- ps_pipiens_whole %>% subset_samples(qPCR_Ratio_Target_Ref >0)

adonis_whole_infection2 <- adonis(vegdist(t(otu_table(ps_pipiens_whole_qpcr)), method = "bray") ~qPCR_Ratio_Target_Ref*Infection*Strain,
       data=as(sample_data(ps_pipiens_whole_qpcr), "data.frame"), permutation = 9999)

res_adonis8 <- adonis_whole_infection2$aov.tab
res_adonis8 <- fill_significance(res_adonis8,"Pr(>F)")

writeData(wb, "NMDS", "Effect of Wolbachia infection (HCA groups), qPCR values and Strain in Culex pipiens with qPCR values", startCol = 1, startRow = 63)
writeData(wb, "NMDS", "ADONIS (Analysis of variance using distance matrices):", startCol = 1, startRow = 65)
writeData(wb, "NMDS", res_adonis8, startCol = 1, startRow = 67, rowNames=TRUE)

saveWorkbook(wb, "Supplementary_Table_2.xlsx", overwrite = TRUE)